JOURNAL OF CHEMICAL PHYSICS 


VOLUME 112, NUMBER 3 


15 JANUARY 2000 


Dipolar relaxation in a many-body system of spins of 1/2 

Alexander A. Nevzorov and Jack H. Freed 3 * 

Department of Chemistry and Chemical Biology, Cornell University, Ithaca, New York 14853 
(Received 25 August 1999; accepted 22 October 1999) 

The method utilized in Paper I [J. Chem. Phys. 112, 1413 (2000)] for treating the density matrix 
equation for a two-spin system in the presence of the dipolar interaction that is randomly modulated 
by translational diffusion, is extended to a many-body system of identical spins of 1/2. Generalized 
cumulant expansions are used, which allow one to take full advantage of the statistical independence 
of the motions of spins. In the high-temperature approximation (appropriate for dilute solutions), for 
a single nonselective pulse, the symmetry of the problem allows one to obtain a compact ordered 
binomial expression for the free-induction decay signal that is related to the two-particle solution, 
and it still contains the two spin-isochromat components. The latter are evaluated by solving the 
corresponding stochastic Liouville equation, which allows one to recover in a unified way the two 
limiting cases including Anderson’s result for statistical broadening in a rigid lattice and the 
classical Torrey-Bloembergen-Redfield expression for the motional narrowing, as corrected by 
Hwang and Freed. The line shape expression in the thermodynamic limit, i.e., for large numbers of 
particles in a macroscopic volume, is obtained. It is found that the many-body dipolar line shapes 
are very close to Lorentzians over the entire motional range studied, with the linewidths 
proportional to the spin concentration, as predicted earlier for the limiting cases. Linewidths plotted 
versus the values of the translational diffusion coefficient clearly show the solid-state limit, the 
motional-narrowing limit, and the intermediate region. The method is extended to describe the 
behavior of the many-body system in a solid-echo sequence. This enables one to obtain the 
homogeneous Tf s over the whole range of motions. A minimum in T 2 is found at approximately 
the same value of translational diffusion coefficient as was found for the two-spin case in Paper I. 
© 2000 American Institute of Physics. [S0021-9606(00)02903-2] 


I. INTRODUCTION 

The problem of the description of spectral line shapes in 
the presence of dipolar interaction amongst a collection of 
spins is well known in magnetic resonance spectroscopy. 
Since there is an r -3 dependence of the magnitude of the 
dipolar coupling on r, the vector connecting any two spins, 
the calculation of spectral line shapes in such spin systems is 
essentially a many-body problem. Indeed, if only a selected 
pair of spins were considered, the equilibrium probability 
density for the distribution in r would be inversely propor¬ 
tional to the sample volume V, thus yielding a negligible 
effect of the dipolar interaction on spectral line shapes for 
macroscopic V. Our goal, therefore, is to consider all the 
spins as a single many-body system and to find the thermo¬ 
dynamic limit for the spectral line shape as the, number of 
spins N —while at the same time V —In earlier work 
by Torrey, under conditions of motional narrowing, the ef¬ 
fect of multiple spins was treated by first solving the two- 
spin problem and then introducing in ad hoc fashion the 
concentration factor rendering the many-body interaction as 
a sum of N two-body interactions for the N+ 1th particle. 1-3 
This yielded Lorentzian line shapes with the spin-lattice, 
= and transverse, R 2 =Tf\ relaxation rates which 
depend linearly on concentration. In solids, magnetic reso¬ 
nance line shapes can be calculated using the statistical 
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theory of Anderson, 4,5 which also yields simple Lorentzians 
over a broad frequency range. However, in the latter case the 
spectra are inhomogeneously broadened. 

In this paper we examine more rigorously the many- 
body effects on spectral line shapes and observed time- 
resolved signals. The emphasis is placed on spin relaxation 
by translational diffusion, which modulates the dipolar inter¬ 
actions between spins. Our general method yields a soluble 
stochastic Liouville equation procedure for the many-body 
problem in the thermodynamic limit, which allows one to 
calculate spectra over the whole motional regime and to re¬ 
cover both the solid-state limit and the motional narrowing 
regime as limiting cases. Our procedure may be summarized 
as follows. In Sec. II, we expand the density matrix in a set 
of eigenoperators of the unperturbed Hamiltonian (cf. Paper 
I, i.e., Ref 6), which enables us to conveniently separate the 
quantum spin variables from the classical stochastic vari¬ 
ables, i.e., Yjj ij = 1JV in the A-spin case. Then we obtain 
a system of coupled differential equations for the expansion 
coefficients, the solution of which can be written in terms of 
an ordered exponential. In Sec. Ill, we show that the 
ensemble-average of the latter can be written in terms of 
connected averages, or cumulants, 7-9 which allows one to 
find the solution for the observed magnetic resonance signal 
in terms of a time-ordered binomial sum of the two spin 
isochromats related to the solution of the two-body problem. 
Finally, in Sec. IV the Markov method 10 is used to find the 
thermodynamic limit for the line shape expression when N 
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—►<*>, F—This leads to an integral equation for the many- 
body signal in terms of the two-spin signal, that was calcu¬ 
lated in Paper I. 6 The integral equation is then solved by 
straightforward numerical methods, (cf. Sec. IV), and the 
continuous-wave (cw) line shapes are described in Sec. VII. 
In this work we use the high-temperature approximation to 
lowest order and do not consider any higher-order effects of 
the sort studied elsewhere. 11,12 This may also be taken as 
implying dilute solutions. Here we show that the cw line 
shape is a Lorentzian over the whole motional range. We 
also show rigorously that the linewidth is linearly propor¬ 
tional to the concentration of spins in the sample over the 
whole motional range. 

The methodology is generalized to the solid-echo experi¬ 
ment in Secs. V and VI. In Sec. VII we investigate how 
motion quenches the echo formation. By simulating the 
SECSY (spin-echo correlation spectroscopy) spectra based 
on the solid echo, we are able to decompose the total Lorent¬ 
zian linewidths into their homogeneous and inhomogeneous 
components. In addition, we show how increasing concentra¬ 
tion quenches the echo even when there is almost no appre¬ 
ciable motion because of the fractional power-law depen¬ 
dence of the homogeneous T 2 on the motional rate. This 
effect is distinct from instantaneous diffusion 4 which reflects 
the fact that a standard i t pulse will not refocus dipolar in¬ 
teractions between like spins. 


II. COUPLING OF THE FREE-INDUCTION DECAY 
COMPONENTS WITHIN THE EIGENOPERATOR BASIS 


Here we consider a system of N particles of spins 1/2 
which diffuse in a sample of volume V. The Hamiltonian for 
such a system in the rotating frame can be written as 

N 

= H {iJ) - (2.1) 

i<j 

Here the double sum is taken over all N spins i and j forming 
distinct pairs such that i<j. For like spins the unperturbed 
Zeeman Hamiltonian H 0 R is given by 

N N 

= 2 (fi - o, rf )4 ;) = AOS /f. (2.2) 

7=1 7=1 


Here we consider the high-field approximation and retain 
only the first term of the dipolar Hamiltonian, namely 


H (U) = 

n R 






[4°/^ - 

(2.3) 


where X =y/(l 6Tr/5)y 2 A is the coupling constant, r t j is the 
distance between the zth and yth spins, H /y - is the orientation 
of the vector r /y - connecting the two spins with respect to the 
main magnetic field. If r /y - and/or Qjj are random functions 
of time, the dipolar interaction gives rise to line broadening. 
The equation of motion for the many-body density matrix is 

^ = -/[// 0 ,p(7)]-72 (2.4) 

at /< j 


It should be noted that the the density matrix contains both 
explicit and implicit time dependence, the latter arising from 
stochastic modulation through the vectors r^. We seek the 
solution of the density matrix equation in the form 

= (2-5) 

M 

where the coefficients g{ e }(0 are functions of time, and the 
are a set of basis operators with numbering scheme (e), 
that are introduced as follows. For a system of N particles of 
spin 1/2, operators can be written as a direct product of 
N spin operators corresponding to different particles 13 
N ' 

£ {f} =n 1*4- ( 2 . 6 ) 

Here {e} denotes a permutation set which provides a certain 
numbering prescription for the eigenoperators, {e} 
= {e x ,e 2 ,...e N }, €j=a, /3, +, or —. I + and/_ are conven¬ 
tional raising and lowering spin operators and the polariza¬ 
tion operators I a and Ip are defined as Ref. 14, cf. Paper I 

4=1+4 (2-7) 

As follows from the trace properties of /+,/_, I a , and Ip , 
the operators of Eq. (2.6) form a complete orthonormal basis 
set in operator space with the Frobenius trace metric defined 
as 

0E {e] ,E [e , ] )=Tr(El ] E {e , } )=S {eh{e , ] , (2.8) 

where the dagger f denotes the Hermitian conjugate. By us¬ 
ing the above orthogonality property of the basis operators, 
one obtains a system of equations which couple the FID (free 
induction decay) components g{ 6 }(/). That is 

^ = — /Afl{ e }g{ e }(/) — C{^{ e ,}F( r / y)g{ 6 /}(z L ). 

(2.9) 

The matrix of frequency offsets is given by [H 0 ,E^] 
= , and the matrix C® is defined by 

^cjJ {e , } F(r iy ) = Tr(£[ e} [// ( 6'U {e , } ]). (2.10) 

The explicit form of clearly depends on the commuta¬ 
tion relations between the basis operators and the Hamil¬ 
tonian describing the interaction among the spin particles. In 
the case of the dipolar interaction in the high-field approxi¬ 
mation, Eq. (2.3), one- obtains the following commutation 
properties for different types of pairwise combinations 

€ i € j 

= [H UJ \l ( i ) I ( i ) ] = 0, 

[i/ ( 67,/W/0)]= +^ F ( ry )(I/(0/0) + (2.11) 

[//< A/W/W] = X F(r tJ )(yVW ~ i/O'/W), 

[//<7),/<;>/67] = X F{ r tJ ) (- }/^/^), 
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where F{Y i j) = Y^{£l i j)lr]j. Thus in the case of the Hamil¬ 
tonian given by Eq. (2.3), the matrices are real and 

symmetric. Out of the possible 4^ time-dependent coeffi¬ 
cients g{ e }(r) of Eq. (2.5), most do not contribute to the 
observable signal. Indeed, the observable free induction de¬ 
cay (FID) signal can be written as 


G(t) = Tr( pg & ] =2 Tr(^ {e} / ( i } ). (2.12) 


Therefore, a nonzero contribution to the observed signal 
arises from the coefficients g{ 6 }(0 corresponding to the op¬ 
erators E{ e } that contain only one operator /+ and. a mixture 
of polarization operators I a and Ip . Others will yield a zero 
trace upon multiplication by an . Moreover, as can be 
seen from the commutation properties, they will be coupled 
via to those coefficients which correspond to operators 
having the total number of operators 7 + that exceeds the 
number of operators /_ by +1. The difference between the 
number of operators I + and /_ in an eigenoperator E we 
shall call the coherence index jul , [H 0 ,E^] = {lAklEfy . 
Thus, in calculating a pure FID signal only the components 
corresponding to jul= +1 are important, (the counter-rotating 
component corresponding to jx— — 1 is related to it by com¬ 
plex conjugation, cf. Paper I). We, therefore, rearrange the 
set { 5 } into four subsets, and we introduce the following no¬ 
menclature for the operators E | e | with jx=l that are neces¬ 
sary to calculate the observed FID signal, Eq. (2.12), from 
Eq. (2.9): 


E 


{*} 


=E 


(Km) 


ri(k,m) 

“maV'W'hii"’}-} 


jd\}. .jdm)jdi\. 1 .. 

a a a ft ft ft + + 




(2.13) 


where m + m 1 + 2k+ 1 = N, m denotes the number of opera¬ 
tors I a , k is the number of pairs /+/_ , set {i} a denotes all 
the possible choices of m particles out of a total N for m 
operators I a , and so on. Two indices k and m are sufficient 
for the classification of the basis operators E and for the 
establishment of the general character of couplings amongst 
the corresponding coefficients g|^j w) (/). Strictly speaking, 
one should also have an additional index jjl corresponding to 
the coherence index of the eigenoperator. However, since 
only the eigenoperators having /z= +1 are important for cal¬ 
culating the FID signal, this index is omitted here. The num¬ 
ber of such operators with given k and m is 


N\ 


N\ 


(2.14) 


NJNp\N+\NJ. m\(N-2k-m-\)\{k+\)\k\ ' 

And the total number of operators having /ul= +1 for 

N particles of spin 1/2 is given by 


2 n ~^ 


N \{2k+^_ u 


2 

k= 0 \ 


2 k + fi)\ k 


2 N \ 
NEjxj 


/*=! 


2 N \ 
N+l)’ 


(2.15) 


which can be proved by induction with respect to N. In terms 
of the above nomenclature the expression for the FID, Eq. 
(2.12), becomes 

N- 1 

G(0=E SgirtO. (2.16) 

m = 0 {*} 

As one can see from the commutation relations, Eq. (2.11), 
the mixture of a, /3, +, and — terms gives rise to the cou¬ 
pling amongst the coefficients g{ k e } m \t) and 
that is amongst those that have the same value of l=k+m. 
Therefore, one can rewrite Eq. (2.9) for the separate FID 
components. gfy m \t) as the following block-matrix diagram: 


(k,m) 



which describes the time evolution of the overall vector of 
coefficients g (t). In Eq. (2.17) we have rearranged the total 
vector g(t) into N parts or sub-vectors, g/(/), /=0,1,...A 
— 1, that decouple from each other, so we can consider them 
separately. Note that only the components that correspond to 
the shaded areas (i.e., with k — 0) contribute to the observ¬ 
able FID signal. 


III. ENSEMBLE-AVERAGED SOLUTION USING THE 
GENERALIZED CUMULANTS 


For the case of like spins, AD| e } = AH for all {e}, one 
can eliminate the first term on the right-hand side of Eq. (2.9) 
by performing the transformation into the interaction repre¬ 
sentation 


(3.1) 

Below we shall solve for the but we shall drop the 

carat symbol for simplicity in notation. One is in fact inter¬ 
ested in an ensemble-averaged solution for the overall vector 
of coefficients, g(r), which is given by an ordered matrix 
exponential for each of its decoupled sub vectors g/(t) 


g/(t) = ( exp o 


N 

~‘x 2 c 

'i</i 


Jo 


F (r ;,/,(M) 


g/(0). 

(3.2) 
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Here the time ordering O is required due to noncommutativ¬ 
ity of the matrices C^; g/(0) is the initial magnetization 
immediately after the first ( tt/2) x pulse, and the sum is taken 
over spins i j and j ] such that i\<j \. In the high-temperature 
approximation, immediately after a nonselective (ttI2) x 
pulse the first nonconstant term of the density matrix corre¬ 
sponding to jjl— 1 is given by 


p,(0) = -2 


A n 

-N _ V j(i)— 9 - 

kT h J+ 




{*} 




(3.3) 


where q = hkl/kT. Thus, for the /th block of Diagram 1 the 
components of the initial magnetization vector g/(0) are 
equal to 1 for and are zero for all others with k 

^0. The number of nonzero components of g/(/) for given 
l — m is determined from Eq. (2.14) and is equal to 


N\ 

m\(N—2k — m— \ )\(k+\)\k\ 


= N 


N- 1 
m 



(3.4) 


The expression for the ensemble-averaged time-dependent 
magnetization vector can be further rewritten in terms of the 
connected averages or generalized cumulants 7 

{oo N N N 

2 (-**)" 2 2 2 

H=1 i \ < j\ i 2 < h ’n^'n 


XC 


(hJ l) 
/ 


^7 


r t 

[>„-* 

dt\ dt 2 " 

dt„ 

Jo Jo 

Jo 


X(f(r Vl (/i))F(r, l/l (^))-f(r fi/i (0)> c ]B(l». 

(3.5) 

where subscript c stands for the cumulant (connected) aver¬ 
age. Here we follow Kubo’s prescription for constructing 
generalized cumulants, in which the ordering O is preserved 
for both cumulants and the exponent to account for noncom¬ 
mutativity of the C-matrices. 

Up to this point our many-body formulation has been 
relatively general and not yet in a form that can lead to 
useful solutions. We now introduce the important assump¬ 
tion that the relative motions of the different pairs of spins 
are stochastically independent, i.e., (Ffoy^fi)) 
^F( r i 2 j 2 ( t 2 ))) = {F( r i ] j ] ( t \)))(k'( r i 2 j 1 ( t 2 ))) unless i, = i 2 
and j j =j 2 . Also we assume that at any time t the spins are 
randomly distributed. This usually means that 
(F(r /i y i (/ 1 ))) = 0 or can be redefined to obtain this result. 
Thus, the cumulants in Eq. (3.5) will vanish 7 unless i\ — i 2 
= "'in and J i =j2 = ' • 'in leading to 

{ °° w 

2 (-*)»2 [C$' J) ]” | dt } j dt 2 ■■■ 

n= 1 i<j JO JO 


In Eq. (3.6) the indexing of pairs has been omitted in the 
cumulants since all pairs behave equivalently. The cumulant 
average over all pairs is thus replaced by an average over one 
pair since motions of different pairs are assumed to be sto¬ 
chastically independent. Actually, to obtain Eq. (3.6) from 
Eq. (3.5) we only require that at least one r f j is uncorre¬ 
lated from the others in each cumulant average. When we 
pass to the thermodynamic limit (i.e., N-+°o) in the next 
section, then even if there is a finite number of correlated 
r f j , their contribution would be expected to be negligible. 

The in Eq. (3.6) still retain the many-body aspects. 
But, it is shown in Appendix A that g ; (0) is an eigenvector 
of the matrix sums for each block in Eq. (2.17) 

with the corresponding eigenvalue given by 

N 

2 [cS' v) ]”g,(0) = [/(|)” + (^- 1 —/)( —i)”] g/ (0). (3.7) 

i<j 

The expression for the FID, Eq. (2.16), can then be calcu¬ 
lated yielding 

G(t)= 2 2 g$\t) 

i-o {*} w 


= 2 N q exp( — iAQt) 



X[iay+(N-l-l)(-±) n ]K n (t) . (3-8) 


(See also Sec. VI where more details are given for the more 
general case.) 

Here we have used Eq. (3.4) to determine the lengths of 
the nonzero components having k = 0 and m = l, 

that yield the relative intensities for each g/(/). Also, the 
ordered cumulant function K n (t) is defined by 

ft f 1 1 [**-1 

K n (t)= dt \ dt 2 dt n 

Jo Jo Jo 

x(F(r(6))F(r(/ 2 ))-T(r(/„))) c . (3.9) 

One can express the ordered exponential in Eq. (3.8) as an 
ordered product of two exponentials, viz. 


expo 


2 (- ix)V (!)"■+ (N- 1 -/)(■- \) n ]K n {t) 


= O exp o\ 


Xexp 0 


1 2 (-ih) n K„(t) 


n- 1 


(N— 1 —/) 2 ( + ih) n K„(t) 

n= 1 


(3.10) 


fj dt n (F(r(t i ))F(r(t 2 ))-■ '7 r (r(/„))) c J g/(0). 

(3.6) 


where the nature of the time-ordering prescription of the 
product may be obtained from expanding and comparing 
terms on both sides of Eq. (3.10) (cf. Ref. 7). It is then easy 
to see that Eq. (3.8) can be rewritten in a binomial form 
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G(t) = 2 N qN Gxp( — /‘Afl/)0j ex Po 


2 {i\xYK„{t) 


+ exp 0 


2 (-ih)"K„(t) 


n= 1 


N- 1 


(3.11) 


Equation (3.11) is a remarkable result, which means that the 
FID for N particles of spins 1/2 can be simply written as a 
time-ordered binomial sum of (N— 1) th degree of the corre¬ 
sponding ensemble-averaged spin isochromats from a single 
pair of particles [compare Eq. (3.11) to Eq. (2.14) of Paper 
I]. The only assumption that has been made is that the mo¬ 
tions of all particles are stochastically independent. In fact, in 
deriving Eq. (3.11) we did not even assume any specific 
form for the F(r(7)) or the detailed nature of the stochastic 
process responsible for its modulation. In treating the 
quantum-mechanical spin variables the only requirements 
were: (i) Identical spins of 1/2 and (ii) a purely secular form 
for the spin Hamiltonian of Eq. (2.1). The particular form of 
the spin-part of Eq. (2.3), relevant for the dipolar interac¬ 
tions, merely determined the eigenvalue of 3/4 that shows up 
in Eq. (3.11). 

To recover the motional narrowing regime as a limiting 
case, one can truncate the cumulant expansion of the expo¬ 
nential in Eq. (3.11) at the second order. Assuming for sim¬ 
plicity that the system is on resonance, AX1=0, and that 
(F(r(0)) —0, this yields the well-known result 2,15 


qN 

G(t)= — O exp£ 


16 


qN 


~ (N— 1) 

X 2 [dt x dt 2 (F{r{t x ))F{r{t 2 ))) c 
Jo Jo 


= T eXp 


9 r t->°° 

dr(t-T ) 

16 Jo 


X(F(r(0))F(r(r))) 


qN 

= T eXP ' 




(3.12) 


where 7(0) = 2/q 7r(F(r(0))F(r(r))) is the zero-frequency 
spectral density for translational diffusion. The explicit ana¬ 
lytical expression for 7(0) has been calculated by Torrey, 1 
and corrected later by Hwang and Freed 3 who have treated 
the excluded volume more rigorously. For the translational 
diffusion with coefficient D T and the distance of minimum 
approach 7, the simplified treatment yields: 7(0) = (2/15) 
X(l/F)(2 /dD T ), cf. Eq. (114'), p. 302 of Ref. 2. In the more 
rigorous treatment, the coefficient 2/15 should be replaced 3 
by 4/27. Here we have explicitly included the sample vol¬ 
ume, V, which arises from the averaging of the two-particle 
correlation function of Eq. (3.12) with respect to the equilib¬ 
rium probability density. Also, we have defined D T as the 
relative diffusion coefficient. When D T is replaced by 2D,. 
where D is the diffusion coefficient for individual particles, 
then the original formula of Ref. 2 is recovered. This shows 
rigorously that under the conditions of motional narrowing, 


r 


the FID signal from N identical spins of 1/2 is given by a 
single exponential, the decay of which is simply proportional 
to the the concentration of spins, C, in the sample, and the 
intensity is proportional to the number of spins, N. 

Another important special case is when motions become 
very slow. In this case the cumulant expansion of the en¬ 
semble averaged exponentials, Eq. (3.11), can be rewritten as 
ordinary exponentials, oscillating with a constant frequency 
3/4^F(r), which are integrated over r. That is 


oo 


exp o\ 


2 

n = 1 


(±4*) K„(t) 




ex Po| 

d 3 r 

T 


3 [‘ 

±i 4 X ) 0 dt ' 


F(r(t')) 


exp[±i|^F(r)r], 


(3.13) 


The resulting FID signal can then be written as 


G(t) = 2~ N qN 


[g-'(3/4)^(r)/ +e /(3/4)^(r)(] 


(3-14) 


One can evaluate Eq. (3.14) in the thermodynamic limit, N 
- 3-00 and K—by rewriting it as 


qN 

G(t)cc lim—— 


c f 

1 _ 2N 2lr J 


sin 6d 0 



r 2 dr 


} N- 1 

X [i— 6 / ( W (r) '] , (3.15) 

so that the concentration of spins C=N/V is kept constant. 
Calculation of Eq. (3.15) yields a single exponential with a 
decay rate of (2'jT 2 /3\/3)y 2 hC, which constitutes the well- 
known Anderson result. 2,4,16 


IV. THERMODYNAMIC LIMIT OF THE MANY-BODY 
LINE SHAPE AS A/-+ 00 , y-+a> 


Equation (3.11) can be evaluated in the limiting case of 
an infinite number of particles in a manner analogous to the 
Markov method used in the ESR of solids. 17 For simplicity 
we shall assume that the system is on resonance, AO=0, and 
rewrite the cumulant expansions, Eq. (3.11), in terms of 
ensemble-averaged ordered exponentials using the first 
equality in Eq. (3.13) 


G(t) = (2- N qN)0\ exp 0 


-xf'dt’FW)) 

Jo 


+ ^exp 0 


;- x J o dt'F(r(t')) 


(4.1) 
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Letting N= 2 we recover the FID signal from a single pair of 
particles which consists of two-spin isochromats, cf. 
Eq. (2.14) of Paper I. If not for the ordering O in front 
of the binomial expression, Eq. (4.1), the functions 
(exp^f+zl^/o^/'F^/'))]) could be immediately evaluated 
by solving the corresponding stochastic Liouville equation as 
described in Paper I. However, for the case of a very large 
number of particles the obstacle of having this additional 
ordering can be eliminated. For this purpose, we shall find 
the statistical limit for Eq. (4.1) at N -*°° while at the same 
time F—As will shortly be seen, the limit F-><*> can be 
treated by introducing a slightly different notation for an 
ensemble average. 

The ordered exponentials are averaged by means of joint 
probability densities. For example, for a stationary Markov 
process, the joint probability density is given by: 

Pr,(r u t\-,r 2 ,t 2 -,...r„,l„) 

= P(r x ,t x \r 2 ,t 2 )---P(,T„-v,t„-\\r„,t„)p ai (r n ,t n ), 


where the prime indicates ensemble-averaging with respect 
to the unnormalized equilibrium distribution. By substituting 
for the volume, V=N/C , one can take the limit A—>°°, so 
Eq. (4.5) becomes 


G(t)= lim— 0\ 1 - 


2NI(WJC) 


X ll-Qxp 0 -i - x 


+ i\-Qxp 0 i ~xj^dt f F(r(t')) 


If the potential U( r) vanishes at |r|—then the integral I 
of Eq. (4.6) approaches unity, and one obtains as a final 
result, in a manner analogous to the Markov method, 10,17 that 

qN C 

G(t)= — Oexp -- 


where p eq (r n ,t„) is the equilibrium probability density, 
P(r„_j ,t n - x \r n ,t n ) is a conditional probability density (cf. 
Paper I). The Boltzmann equilibrium distribution is given by 


Peq(rT) = Peq(r) 


exp[— U(x,y,z)/kT\ 


/o*/o' v /n Zjx dy dz exp[- U(x,y,z)/kT] 


X / 1 — exp 0 -i - xj Q dt'F(r(t')) 


+ ^1—exp 0 p ~xJ Q dt'F(r(t')) 
where we have used the fact 18 that 


— lim/(/V) 
= gAf—>“ 


Here we are interested in just the stationary /? eq (r) consistent 
with thermodynamic equilibrium for the many-body system. 
By substituting the variables, the normalization coefficient in 
Eq. (4.3) can be recast in the following way: 


provided that the limit of /(A) exists and is finite. 

Equation (4.7) is still written in a rather symbolic man¬ 
ner which contains ensemble-averaged ordered exponentials. 
It can further be rewritten in the form 


dx dy dz exp[— U(x,y,z)/kT] 


= v(' f 1 f 1 d£d v d£exji-U(L x lL y i l ,L I £)/kT] 

Jo Jo Jo 

= VI(L x ,L y ,L z ), (4.4) 

with ij=x/L x , etc. Using the expression for the ensemble- 
averaged time-ordered exponential [Eq. (3.1) of Paper I] one 
can pull out the normalization integral from the ensemble 
averaging in Eq. (4.1) and write that 


G(t) = (2~ N qN)0\ 2- 


VI{L X ,L V ,L Z ) 


X 11 — exp 0 -i - X jyt'F(r(t’)) | 

+ |l—exp 0 d- X ^dt’F{r(t')) J j , 


qN f dg(t') 

qN \ f' dg(t') , . 

= — exp C dt' , , (4.9) 

2 [Jo dt 

where g(t) = [g+(t)+g-(t)]/2, and the individual spin iso¬ 
chromats are given by 


g±(0 = (exp o +i - xj Q dt'F(r(t')) 


= ( exp -hi - xj Q dt'F(r(t f )) 


Correspondence between Eqs. (4.7) and (4.9) can be estab¬ 
lished by differentiating with respect to time the formal se¬ 
ries solution given by Eq. (3.4) of Paper I. By comparing the 
expression for G(t). Eq. (4.1) with its thermodynamic limit , 
Eq. (4.9), we have found that in the limit of a large number 
of particles the time ordering in G(t) becomes unnecessaiy. 

Note that G(t) decays to 0 when t— >co, as. it should, 
since it describes an FID signal. Indeed 
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G(co) = exp 



dg(t') 

dt' 


We can now consider the kernel of the integral equation as 
just a delta function, which yields 


= exp[C(-g(0)-/ftjg(o))| (0=0 )] 

= exp[ —Cg(0)], (4.11) 

where g(<y) is the Fourier transform of g(t ). But g ± (0) is 
equal to the integral of the unnormalized Boltzmann equilib¬ 
rium distribution over volume V, cf. Sec. Ill of Paper I, 
which yields at TV—and F—>°° 

G(°o)= lim exp[ — Cg(0)] 

N -* °° 

V-+CX> 

[ c 

= exp ~2 

= lim exp — C I d 3 rexp[— U(r)/kT]\ 

N —> co l J J 

V-+°° 

= lim exp[ — NI( \[v) ] = 0, (4.12) 

N —* 00 
V—>°° 

where the integral / is defined in Eq. (4.4). Equation (4.9) 
can be rewritten in terms of a differential equation for G(t), 
viz. 


~ / N G(0) 

G(w) ~[Cg(0)/x']-/o,’ 


(4.17) 


i.e., the many-body spectral function is given by a Lorentz- 
ian with half-width Cg(0)/r f , or in the absence of any po¬ 
tential CV/t' =N/t ', where TV—►<» and r'—(or F— 
and r'—>°°), cf. Eq. (4.10). Equation (4.17), of course, de¬ 
scribes the line shape only qualitatively, since one still has to 
solve for r as a function of volume F and the diffusion 
coefficient D T until one converges to the thermodynamic 
limit. 


V. TREATMENT OF MULTIPLE PULSES USING THE 
EIGENOPERATOR BASIS: SOLID-ECHO 
SEQUENCE 

In this section, we consider effects of multiple nonselec- 
tive hard pulses applied to a many-body spin system. As an 
illustrative example, let us consider a solid echo experiment 
involving an intermediate (tt/2) v pulse, viz. 


dt dt 


(4.13) 


By performing the Fourier-Laplace transform of Eq. (4.13), 
one can write that 


~ iC f+°° ~ 

ia)G(a))— -— I (co — o)')g((i) — (L>’)G((jd')d(i)' 

2 77 J — co 

= [Cg(0)-l]G(0). (4.14) 


- \—j -acquire. (5.1) 

x ' ' y 

Immediately after the intermediate pulse, the density matrix 
in the rotating frame becomes 



p(t+) = 2 g {e} (T)R(0)E [e] R-\e), (5.2) 


Thus, it is sufficient to know just the two-spin line shape 
function g(co) and the concentration C to solve for G(co). 
Equation (4.14) can be transformed into a matrix equation by 
discretizing the convolution kernel with respect to co and co ', 
and then inverted numerically. 

It is easy to see that the solution of Eq. (4.14) is expected 
to be close to a Lorentzian. First we note that the two-spin 
spectral function g(co) is expected to be much narrower than 
the TV-body G(co), and therefore, it is not unreasonable to 
approximate it by an infinitesimally narrow Lorentzian, viz. 


where R(6) is the rotation operator. In the rotating frame, 
R(9) is equal to expwith 6=irl2 for the case of 
solid echo. In general, the rotation of an eigenoperator 
will yield a linear combination of the eigenoperators E{ e '} 
that form a complete basis set 


R(e)E [e] R- l (e)=J J X {e][e , } E {e , } . (5.3) 

{«'} 


g((D-w')*>g(0) 


l — i(co— co 


(4-15) 


The elements of the pulse propagator matrix X T in the eigen¬ 
operator space can be defined as 


(Here we are ignoring the infinitesimal “Pake doublet” for 
simplicity.) Substituting Eq. (4.15) into Eq. (4.14) and using 
the identity flZG(co')dco' = 2 ttG(0), we obtain 


icoG(co) — 


Cg{ 0) 

2 ttt' 


1 — i(co— co')t' 


G(co')dco' 


~-G( 0). 


(4.16) 


^ {f } {e '}-Tr[^ { V } 7?(^ {f }7?“ , (^)] 

= Tr[E {£] R-'(8)El, } R(d)]. (5.4) 


To relate the components g{ e }(0 before and immediately af¬ 
ter the pulse,, we write a. new expansion for the density ma¬ 
trix p(r+) in terms of the eigenoperators 
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p(r+)='Z g{ e }(.r+)E [e} . 
ft} 


where the equivalence of Eq. (5.8) to Eq. (5.7) may be 
(5.5) shown by expanding out the exponential operators. 


Comparing Eqs. (5.2) and (5.5) and making use of Eq. (5.3) 
we find the relation between the vector of coefficients 
g(r+) and g(r) 

g(r+) = Xg(r). (5.6) 


After the pulse, the vector of coefficients continues to evolve 
as given by Eq. (2.9) for each element, but with a time shift 
of t — t+ and with the initial conditions at r+ determined by 
Eq. (5.6). Therefore, the overall vector g(/) in the rotating 
frame, at times t> t+ is given by 


g(0 




XXe 


-/AOr 


)g(0). 


(5.7) 


where AH is the offset matrix. Note that in the above ex¬ 
pression the pulse propagator X in general may mix all 
blocks of the C-matrix, including the components with dif¬ 
ferent coherence indeces /jl as opposed to a pure free- 
induction decay. As will shortly be seen, a nil pulse results 
in mixing of the components with jjl= ± 1. Now, the diagonal 
matrix AH containing frequencies /x(H— corf) commutes 
with the block-diagonal matrix C, since AH is just a con¬ 
stant matrix for each /xth block of C, and we can write that 


gW = e -/Afi(,- T ) (e -<</ V| („)) 


Xe 'i 
o 




)Xe~' AilT g(0), 


(5.8) 


J 


VI. MULTIPLE PULSES: TIME ORDERING AND THE 
GENERALIZED CUMULANT EXPANSION 

The next step is to perform the ensemble averaging of 
the ordered exponential of Eq. (5.8) by taking into account 
the statistical independence of the motions of spins, as has 
been done for the case of a pure FID. Using the standard 
procedure for constructing generalized cumulants, 7 one ob¬ 
tains that 


Qxpo 


X expX 


* f' 

-i*2 dt x F(r, • (f,)) 

M</i 

N r r 

-i*2 Xd'V.*- 1 \ T dt x F{r ,,,(/,)) 
h<J\ 

" [i 

iX 2 C ( f ' 7 .) *, F ( r v (/,)) 


= exp 0 ^expc 


+ XC ( 'Vi>X‘ 




('.)) 


1 , ( 6 . 1 ) 


where c stands for the cumulant averaging. Note that the 
time ordering O also takes care of the fact that C (/ J ^ and 
XC^’/^X -1 d 0 no t in general commute. Indeed, integrals of 
the type fldt^j^dt' cannot be time-ordered since in this 
case t'^t" everywhere, and thus C (/ ,y ^ can never be to the 
right of XC (7,,J,,) X~ 1 after the ordering is performed for the 
second time. If the motions of spins are statistically indepen¬ 
dent, by expanding the second ordered exponential in Eq. 
(6.1), one obtains 


N N N 


N 


exp 0 2 (~*)”2 2 •••2 2 C ( ''i7i)...C ( Wp ) XC ( V+i7 P+ i)...C ( 'V n )x- 1 

»= 1 P = 0 h<J\ i p < jp i p + l < Jp+l * n < Jn 


X 


\ T dh "'\r dtp L dfp+1 '"L" ' dt ”( F(r >Pl (tl ^''' F(r yp (f ^ F( %^Jp+> (t P + ' y) ''' F(r >’d t P ) ^ 


= expo2 (-!*)"2 2 [C (y >P’X[C (y) ]" - ^X _1 

n— 1 p = 0 i<j 


J dt,---J P 'dt p j°dt p+ r"J g dt„(F(r(t,))---F(r(i p ))F(r(/ p+1 ))---F(r(t„))) c 


( 6 . 2 ) 


Using the ordering prescription, one can obtain the expansion of the ordered exponential up to various orders m. Restricting 
ourselves for illustrative purposes to m — 0,1,2 one obtains from Eq. (6.2) that 
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exp o 


00 n N C t C t Ct Ct 

■ 2 (-*•*)" 2 2 [C (y) y’X[C (? -' ) ]” _ ' p X -1 dt x - P ~'dt p T dt p+ r- "~'dt„ 

n = 1 p = 0 i<j J t J t JO Jo 

X (F(r(tj))- ■ ■F(r(t p ))F(r(t p+1 ))- ■ ■F(r(t„))) c X exp(-/Aftr)g(0) 

n N 


= i + 


2 (-/*)" 2 2 [c^Fxtc^r-'x- 1 [dt x -~ [ P ~'dt p [ T dt p+l -~ 
n -1 p=o i<j Jt Jt Jo 

00 00 

f'"' 1 * n (F(r(/i))-”^(r(g)F(r(^ +1 ))---F(r(0)> c + 2 2 

Jo w ' = 1 » = 1 

x 2 2 2 [c (, ''- / ' ) ]"'[c (y) ]" _ - p xtc (v) p’x _1 ['<*',••• \ ,,n '-'dt r n , f \'"~ p ~'dt„_ p f T dt „_ p+I ••• 

P~0 j'cj' i<j Jt Jt Jt Jt JO 

oo oo n' N N 

|'"' 1 *„(F(r( ? ' 1 ))-- J F(r(?'„,)))c(^(r(/ 1 ))-- J P(r(0)}c+2 2 2 2 E 

JO ;7 ' = 1 w 1 p' = i /' <y' i<j 


f ” l o'/„( J F(r(?' 1 ))---F(r(?'„/))) c (/ ; ’(r(?i))---F(r(/„))) ( 
Jo 


H-| X exp( — /Aflr)g(O). 


(6.3) 


The consistency of the ordering application can be checked 
by letting X=l, which corresponds to a pure FID signal. In 2 [cJ^^XfC^] 77- ^ exp( —/ADr)g/(0) 
this case the previous time-ordered ensemble-averaged expo¬ 
nential, Eq. (3.6), is recovered. The latter essentially follows 
from applying the following integral identities: 


f<J 


= cos(ADt)[/(|F(-!) w ^ 


J‘dt' + J T dt' = J'dt', 

\'dt' i' dt"+ i'dt' \ T dt"+ \ T dt' (' dt" 

Jt Jt Jt JO JO JO 

= i'dt' ['' dt", 

Jo Jo 

J'dt'J'dt"J '"dt"' + J'dt' J'dt" J 7 dt”' 


(6.4a) 


+ (Ar-l-/)(-fr(f)"^]g,(0) 

= cos(Afir)(!)”[(— \) n ~ p l+{ — l) p (N— 1 -/)]g,(0), 

(6.5a) 


xg(0) = g(0). 


(6.5b) 


(6 4b) Validity of Eq. (6.5b) also follows from the fact that a pure 
oscillating signal, which has the initial magnetization along 
the y axis, is not affected by an intermediate pulse around the 
y axis in the rotating frame. We shall further consider the 
case when Afl=0, i.e., when the system is exactly on reso¬ 
nance, and define the following integral operators: 


+ J'dt' I 

M 

'"dt"'+ \ T dt' i' dt" [‘"dt"' 

f T 

Cl 


'o J 

o Jo Jo Jo 

I\ = {-i\x)\ dt \> 

Jo 

h=(~ijX) J dt,, (6.6a) 

= i'dt' 1 

\‘ dt" 

i'dt'", (6.4c) 



JO J 

0 

1 0 

Cl 

r t ri 




dt 

{-‘aX) dti = ( — ijx) 1 dt x s(t x ). 

and so on. 



J T 

Jo Jo 

In Appendix B it is shown that, for a nonselective hard 


(6.6 b) 


pulse, each subpart of the starting vector-coefficient g/(0) is 

an eigenvector of the matrix sum 2 / <y[C ( V)]/ 7 x[C (/, - / )] ?7 P } where the ^-function is defined as follows: $(/) = +1 if t 

where / denotes the number of operators I a in the corre- >t and s(t) = — 1 if 0 <t<r. The following short-hand no- 


sponding eigenoperator subset is . That is to say 


tation is also introduced: 
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{l}"(F{t x )F{t 2 y--F{t n )) c 

■if - '?*) //''/.V 


Xs(t l )s(t 2 )- ■ ■s{t„){F{t x )F{t 2 y ■ ■F(t„)) c 



” f f'l f'n- 1 

dt i | dt 2 ‘ dt„ 

Jo Jo Jo 


Xs(t,)s(t 2 )- ■ ■s(t„)(F(tOF(t 2 )-■ -F(t„)) c , (6.7) 

which can be further rewritten as [compare also with Eqs. 
(6.4a)-(6.4c)] 

{lY(F(t,)F{t 2 )---F{t n )) c 

= (-'lx>''E \dt\---[ P dt f dt +i --- 

p = 0 J T Jr Jo 



dt n {-\y-p{F{t x )F{t 2 y--F{t„)) c 


o{iy-■ ■ 

m 

= oU {/}"'• 

7—1 

" m + “- + n 2 + n x 

= E 0{I 2 }”” + ■ + "2 + "t-P(-l)P{] l }P. (6.9) 

p—0 


We shall now prove by mathematical induction that the mth 
term of the expansion, Eq. (6.3) can be factorized as 

[/ + (- 1 1 - /)]• ••[/+( — 1 )” 2 (N- 1 - /)] 

X [/ + (- 1 )"i (N- 1 - /)]0{/}' 7 - • -{/} 77 2{/} /7 i 

m 

= oU [/ + (-l)"/(A7-l-/)]{7}"/, (6.10) 

7=1 


n 

= o E {hY~ p i - 1 ) p {yy{F{t x )F{t 2 y - ■F{t „)) c . 

P = 0 

( 6 . 8 ) 

The ordered product of integral operators acting on various 
cumulant partitionings n l ,n 2 ,...,n m can also be written sym¬ 
bolically as [cf. Eqs. (5.15a)-(5.15c)] 


from which the ordered cumulant exponential function can 
be recollected later. 

The m = 0 and m = 1 terms are evident from the eigen¬ 
value properties of the starting vector g/(0). Applying the 
ordering for the m+ 1th term of the expansion and operating 
on the starting vector g/(0), and using its eigenvalue prop¬ 
erties, we obtain 


J 


OO 00 00 00 

E E -S E 

n m +1 = 1 n m= X ”2= ] ”1 = 1 


( N m 

[/+(-1 )"-+ >(n- i - i)]o{i 2 y^ in [/+(-i y<(N-1 - /)]{/}"- 


”m +1 "m 

+ E [(-i)./ + (-ir--^-i-/)]o { / 2 }— { /,Fn [(-i)”'7+(t/-i-/)]{/}"' 

p= 1 7=1 


( 6 . 11 ) 


where N m =n x + n 2 -\ -h n m . Substituting p—>p — N m , in the second term and using the notation of Eq. (6.9), we can add up 

the two terms yielding 


00 00 00 00 f m + 1 N m 

E E -E E II [/■+(-! y>{N-\-i)-]o{i 2 y ^E 

”m +1 = 1 n m =) - n 2 = 1 ”1 = 1 { i=l P = ° 

N m +1 m 

+ E [/(—i y~ Nm +( —i ) Nm+] ~ p (N— i—/)]E [/(-1 )"'+(#- \-i)~\o{i 2 } Nm+ '~ p {p} p 

P = N m +\ 7=1 

777 + 1 N m + I 777 +1 

= n [/+(-iw-i-/)] E o{i 2 } N ^-p(-mhY=oU [/+(—i) n '(Af— i—/)]{/}"'■, 

7=1 p — 0 7=1 


( 6 . 12 ) 


which completes the proof by induction. Collecting back the ordered cumulant exponential function, we get 


O exp E [l+(-iy , (N-l-l)]{iy(F(t l )F(t 2 y--F(t„)) c 

n=\ 

00 

= ex Po E (-/+)”[(!)"/+(-fr(A?-l-/)] [dtA''dt 2 -- f''" , r/ V (r I >(7 2 )--^(/J(^(7,)7 ? (/2)--^(^))c- (6-13) 

77=1 JO JO JO 
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Thus, one can rewrite Eq. (5.7) in terms of usual ensemble-averaged exponentials by using Eq. (6.13) 


exp o 


-ixH C ( 'V.T' 

i\<J l Jt 


dt i F ( r i l jS t J) 


XexpQ 


X exp 0 


3 f ' 

~i 4 X] dtMti 


)F(r (?])) exp 0 | 


r 

~iX 2 C('V.) T ^,F(r, ,(/,)) 

fi</i Jo 171 

n-i-i 


g/(0) = O 


' i *fo 


dt,s(t,)F(r ((,)) 


g/(0). 


(6.14) 


Summing over the components having total statistical weight of A , (' v / '), Eq. (3.4), one can finally rewrite the 

observed solid-echo signal in the rotating frame on resonance 


G(t) = (2~ N qN)0 exp 0 


JX f dt's(t')F(r(t')) + exp 0 i jxf dt' s{t')F{r{t')) 
4 Jo 4 Jo 


N- 1 


(6.15) 


i.e., a binomial expression is again obtained, as in the case of 
a pure many-body FID signal, Eq. (4.1). The latter is easily 
obtained from Eq. (6.15) by setting ,?(/)= 1 everywhere. 
However, the derivation of Eq. (6.15) could be also regarded 
as a formal check of the consistency of the introduced time¬ 
ordering “O” for treating the ensemble averaging of expo¬ 
nential operators, which in general do not commute but are 
applied during different periods of time, i.e., one is from 0 to 
r and the other is from r to t. 

Equation (6.15) can be extended to the limit of a large 
number of particles exactly in the same way as it was done 
for the FID signal by using the Markov method. Taking the 
limit F—>°° yields 


qN 

G(t)= — exp' 


C 

2 


1 -exp 0 


+ 1 -exp 0 


~i ^ X J o dt's(t')F(r(t’)) 

i ^ X J o dt's(t')F(r(t')) 


(6.16) 


which is to be compared with Eq. (4.7). The ensemble- 
averaged exponentials are evaluated by solving the appropri¬ 
ate modified stochastic Liouville equation taking into ac¬ 
count the ^-function as described in Paper I [cf. Eqs. (4.9) 
and (4.10) of Paper I]. 

For very slow motions, i.e., as Z> 7 —>0, in averaging the 
exponentials in Eqs. (6.15) and (6.16) we can hope to ap¬ 
proximate the two-spin stochastic Liouville operator as L ± 
— l/T ,/ Tzf / Y J F(r), where r is an effective diffusional re¬ 
laxation time with (cf. Fig. 5 of Paper I where 

#=1/2). (However, more generally we can expect t to de¬ 
pend on the value of rand r, cf. Ref. 19 for related cases.) In 
this limit, the homogeneous line broadening is expected to be 
much less than the inhomogeneous line broadening. The 
solid-echo signal could then be written approximately as 


qN I N \ I \t-2r\\ 
G(0 a ^-exp|-^7/jexp|-7^“]' ^ 6A7 ' > 

Here the inhomogeneous T *~ 1 =R* is given in the thermo- 
dynamic limit by the Anderson formula, cf. Sec. Ill, and the 


ratio N/t" may be regarded as a homogeneous relaxation 
rate, T^ X in the limit TV—and r"—I n Eq. (6.17), the 
inhomogeneous line broadening gets refocused, whereas the 
homogeneous line broadening yields the decay of the echo 
amplitude. From the behavior of the echo signal described 



- 10 'VVtP o 10 'Vvd 3 

FIG. 1. Many-body line shapes for spins of 1/2 in the thermodynamic limit 
as a function of the diffusion coefficient, D T (in units of y 2 fild), and con¬ 
centration of spins, C of 1.2, 2.4, and 4 (in units of 10 16 cm -3 ) from top 
down. The values for C and D T are as indicated. The ratio r max Id was set to 
100, after checking that this led to convergent results; then C was varied by 
varying the number of spins, N. Concentrations were calculated from C 
— NIV , where F= (47r/3)r^ ax , where d was arbitrarily set to 10 A. Line 
shapes are close to Lorentzians over the entire motional range for all con¬ 
centrations considered. Increasing the concentration of spins yields a pro¬ 
portional increase in spectral linewidths; whereas increasing the diffusion 
rate D T results in motional narrowing. 
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FIG. 2. Line shapes and spectral linewidths, R * (inset) in the thermodynamic limit as a function of the diffusion coeffcient D T calculated for the concentration 
of spins C=2.4X 10 16 cm -3 with d=\0k. Line shapes have been fit to a Lorentzian function, from which the line broadening term, R* , has been 
determined. The plot of R * vs D T (inset) shows the behavior of linewidths over the whole motional range considered. 


qualitatively by Eq. (6.17), one can see that even in the limit 
of very slow motions, one can expect a decay of the echo 
amplitude in a multiple-spin system, i.e., the many-body 
echo may not be refocused completely at 2r. However, be¬ 
low we shall only use the more general expression, Eq. 
(6.16) without assuming any simplified limiting form of the 
type of Eq. (6.17). We shall indeed find that the above 
simple argument has some qualitative validity, but more 
quantitatively 1 is given by a fractional power law some¬ 
what different from D l j 2 . 

VII. SIMULATIONS OF MANY-BODY SPECTRAL LINE 
SHAPES AND ECHOES IN THE PRESENCE OF 
TRANSLATIONAL DIFFUSION 

Figure 1 shows the line shapes in the thermodynamic 
limit associated with a large number of particles allowed to 
diffuse in a large volume. The ensemble-averaged ordered 
exponentials have been calculated by the stochastic Liouville 
equation method (cf. Paper I). Convergence of the spectra 
has been checked by increasing the number of particles, N 
and the ratio r max Id, while keeping the concentration, C con¬ 
stant. [Here we may use dimensionless units for C, i.e., C 
= N/V, where V= (47rl3)(r max !d) 3 ]. It was found that if the 
number of particles N is more than 100 T then an r max /d of 
about 100 is sufficient to obtain a converged line shape. One 


also needs to discretize the convolution kernel of Eq. (4.14) 
with 4000 frequency points. Interestingly, the r max /d re¬ 
quired for convergence was not found to depend appreciably 
on D t , implying that the significant volume of spins sur¬ 
rounding a test spin remains the same, independent of mo¬ 
tional rate. Line shapes are very close to Lorentzians over the 
motional rates considered. This was checked by fitting them 
to Lorentzian functions (see below). Increasing the concen¬ 
tration yields a proportional increase in widths. 

To further explore the many-body line shape function, 
spectral linewidths were calculated as a function of the dif¬ 
fusion coefficient at a fixed concentration, and are shown in 
Fig. 2 (C=2.4X 10 16 cm -3 , J=10A). The linewidths were 
detennined by nonlinear least-square fits of the line shapes to 
a Lorentzian function. As can be seen from Fig. 3 plotted on 
a semilogarithmic scale (inset), the curve has the form of a 
sigmoidal function with a characteristic point of transition to 
the “quasi-solid state” limit at a value of the diffusion co¬ 
efficient D T of —0.05 (in units of y 2 hld). For electron spin¬ 
bearing molecules for which the distance of closest approach 
between their centers is d— 10 A, this value of D T converts 
to 1.64X 10~ 7 cm 2 s -1 , which is typical for spin labels in 
moderately viscous liquids (e.g., liquid crystals 20 ). For C 
= 2 AX 10. 16 cm" 3 one has^a value of R% of 2.61X10 4 s _ 1 . 
One would need higher concentrations to obtain an observ- 
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time, t 

time at which pulse is applied is 10 times shorter 


FIG. 3. Simulation of the solid echoes for N spins of 1/2 in the thermody¬ 
namic limit as a function of diffusion coeffcient D T (in units of ^h/d) and 
concentration of spins, C, of 2.4, 4.8, 12, 24 (in units of 10 16 cm” 3 ) from top 
down, at exact resonance. As C increases, the signal decays faster, yielding 
no echo at higher concentrations. However, by decreasing the time t, at 
which the second pulse is applied, one recovers the echo again (cf. inset). 
The abscissa is in units of time from 0 to 2X10 5 d 3 /y 2 h; inset from 0 to 
2X10 4 d 3 /y 2 h. 


able effect. Thus, increasing C by a factor of 200 (8 
X 10“ 3 M) yields an R 2 °f 5.22 X 10 6 s _1 , which should be 
measurable. 21 In the case of proton NMR (nuclear magnetic 
resonance), D T of 0.05 translates into 3.78X 10 _13 cm 2 s _1 . 
This would mean that slow motional effects set in only for 
extremely slow motions characteristic of high-molecular 
weight polymers or diffusion in solids. 

Figure 3 shows many-body solid echoes at exact reso¬ 
nance. The distance of minimum approach has been taken to 
be d= 10 A. For illustrative purposes only, we have included 
a nonzero carrier frequency. Increasing the concentration of 
spins yields a decrease in the echo amplitude. Thus, going 
from C=2.4X 10 16 cm -3 to C=2.4X 10 17 cm -3 results in 
almost no detectable echo signal in the latter case even in the 
very slow motional regime (Z) r = 1CT 4 ). Increasing the dif¬ 
fusion rate also yields a decrease in the echo amplitude, so 
that it becomes impossible to refocus the dipolar interaction 
for faster motions. This happens for smaller D T as C gets 
larger. Refocusing of the echo at larger C can still be 
achieved by applying the second (tt/2) v pulse much earlier. 
For instance, we find that in order to get an appreciable re¬ 
focusing of the echo signal near the rigid limit, for a concen¬ 


tration of electron spins of C=5Xl0 18 cm~ 3 (or ~8 
X 10” 3 M), the second pulse must be applied no later than 
600 ns after the first pulse (not shown), which is easily 
achieved by pulsed ESR methods. However, this should be 
regarded only as a rough estimate, since in the present analy¬ 
sis additional factors leading to a loss of refocusing (e.g., 
additional homogeneous line broadening, noise, etc.), have 
not been taken into account. 

The above effect of concentration on the refocusing of 
the echo is illustrated in Fig. 4. Here the echo envelopes, i.e., 
intensities measured vs t x = 2r , are plotted at various con¬ 
centrations ranging from C=2.4X 10 16 cm -3 to 2.4 
X 10 18 cm -3 (d= 10 A). As can be seen, at high C the decay 
of the echo envelope is much faster, thus resulting in no echo 
at sufficiently long r s. (The actual procedure we used was to 
calculate the SECSY signal, then Fourier transform it with 
respect to t 2 , and then to study the dependence of the 
ol >2 = 0 slice, cf. Paper I, Sec. V.) 

Figure 5 shows the homogeneous relaxation times T 2 as 
a function of the diffusion coefficient D T plotted for different 
concentrations of spins C. The T 2 values have been obtained 
by fitting the echo envelopes (cf. Fig. 4) to single exponen¬ 
tial functions. This method of determining T 2 s corresponds 
to the SECSY technique as described in Paper I [cf. Fig. 4(a) 
and its discussion in that paper]. The data in Fig. 5 have been 
rescaled with respect to concentration, so that a departure 
from linear dependence on concentration (of about 10%) can 
be seen. This occurs in the slow motional regime, which also 
shows a limiting dependence on the diffusion rate that goes 
as Dj 032 (vs the D ^°' 5 dependence found for a two-spin 
case in Fig. 5 of Paper I). By contrast, in the motional nar¬ 
rowing regime homogeneous relaxation times depend lin¬ 
early on both the concentration, C and D T . Note that in¬ 
creasing C does not shift the T 2 minimum. In addition, the 
T 2 minimum for a system of multiple spins is shallower than 
that for a two-spin system (cf. Fig. 5 of Paper I). 

Finally, a comparison of the spectral linewidths, R *, 
with the decay rate of the echo envelope, R 2 , is shown in 
Fig. 6 (C=2.4X 10 16 cm -3 , d= 10A). Here the log-log 
scale has been used to show the solid-state limit, the inter¬ 
mediate region, and the motional-narrowing limit. In the fast 
motional limit (larger D T values), line broadening is fully 
homogeneous^ and R 2 approaches R* . The curve then be¬ 
comes linear with a slope of — 1 which corresponds to the 
motional-narrowing, or Redfield limit as given by Eq. (3.12). 
The intermediate motional limit occurs near the position of 
the R 2 minimum, i.e., ~D T = 0.05, and spans about two or¬ 
ders of magnitude in D T . As D T — >0, R 2 values approach 
the limiting Anderson expression for linewidths given by Eq. 
(3.15), whereas R 2 further decreases with D T . Thus, in the 
slow motional limit there is a clear separation between R 2 , 
which mostly consists of inhomogeneous line broadening, 
and R 2 , which may be regarded as homogeneous line broad¬ 
ening. Note that at D T = 0.05, R 2 ^R* ; thus, for the above 
example of electron spins with C^8X10~ 3 M, R 2 would be 
about 5 X 10 6 s" 1 , which can be a significant contribution to 
the homogeneous linewidth. 19 
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FIG. 4. Effect of concentration on echo refocusing. Here the echo envelopes are calculated at various concentrations as indicated for a very slow motional rate 
( D t = 10 ~ A y 2 fild). The intensities are measured from a SECSY simulation [cf. Paper I, Fig. 4(a)] that is Fourier transformed just with respect to t 2 . Then 
the slices at £02 = 0 were plotted as a function of the delay time =2r. Increasing concentration results in no echo at sufficiently long r s. 


VIII. SUMMARY AND CONCLUSIONS 

The main result in this work is the solution for the mag¬ 
netic resonance line shape arising from N identical particles 
that bear spins of 1/2 and are undergoing dipolar interactions 
which are modulated by their relative translational motion. 
This result describes the whole motional range from the rigid 
limit through the motional narrowing regime. The general 
result for the FID, Eq. (4.1), only required the assumption of 
the stochastic independence of the motions of the spins and 
the usual high temperature approximation, Eq. (3.3). Equa¬ 
tion (4.1) shows that the many-body solution for N identical 
spins of 1/2 may be obtained from just a knowledge of the 
behavior of a single pair of particles, and it is written as a 
time-ordered binomial expression involving the two-spin iso- 
chromats associated with the two allowed transitions result¬ 
ing from the interacting pair of spins of 1/2. In the very slow 
motional limit, the line shape expression yields Anderson’s 
statistical theory for line shapes in a rigid lattice. In the mo¬ 
tional narrowing regime it is consistent with the well-known 
Torrey results. In the thermodynamic limit, i.e., when N 


—and F—>°°, under the assumption of a stationary process 
(i.e., a system at equilibrium), a simple integral equation, Eq. 
(4.14), has been obtained for the many-body line shape func¬ 
tion, G(co) in terms of the two-spin line shape function g(co) 
and the concentration, C=N/V. 

The general expression for the many-body FID signal, 
Eq. (3.11) was generalized to the proto-typical solid echo 
pulse sequence, Eq. (6.15), and then to the thermodynamic 
limit, Eq. (6.16). The solution was rigorously developed for 
the case of exact resonance. This case guarantees that no new 
orders of coherence besides ± 1 are produced by nonselective 
pulses in a system of N identical spins of 1/2. (This follows 
rigorously from the eigenvalue properties of the starting vec¬ 
tor and the interaction matrices in the eigenoperator repre¬ 
sentation.) This result also depends on the high-temperature 
approximation for the equilibrium density matrix, as well as 
the assumption of the statistical independence of the motions 
of the spins. Removing either of these assumptions would 
result in_ generating other orders of coherence,., as is well- 
known in solid-state NMR. 22 
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FIG. 5. Homogeneous relaxation times, T 2 , for//spins of 1/2 in the thermodynamic limit as a function of the diffusion coeffiecient D T , calculated for three 
different concentrations as indicated. The values at the higher concentration were rescaled accordingly in order to illustrate the departure from the linear 
dependence of T 2 on D T , as opposed to T* . A T 2 minimum is seen at ~ D T = 0.05 y 2 hld, i.e., at about the same value as for the case of two spins of 1/2, 
cf. Fig. 4 of Paper I. However, for multiple spins, the T 2 minimum is shallower. Dashed lines show fits to limiting power-law behavior of T 2 on D T in the 
slow and fast motional limits. 


Convergent calculations of the many-body line shapes 
and echoes could readily be performed based on the stochas¬ 
tic Liouville equation for two spin-bearing particles under¬ 
going simple Brownian relative translational diffusion. These 
calculations showed that a Lorentzian line, whose width is 
proportional to C, is obtained over the whole motional range 
and is in agreement with the Anderson solid-state limit and 
the Redfield motional narrowing limit. The solid echo simu¬ 
lations show characteristic echoes for very slow motions, 
which are attenuated as the motional rate increases and/or C 
becomes larger. In particular, the homogeneous r 2 ’s, which 
are associated with the echo decay show a power-law depen¬ 
dence in the very slow motional regime (as D ^°' 32 ). Al¬ 
though some departure from a linear dependence of 1 on 
C is found at slow motions (unlike the case for T* _1 ), simi¬ 
lar values are found for the D T ’s that correspond to the T 2 
minima. 

The present work may thus be seen as a unifying method 
which bridges limiting theories of many-body magnetic reso¬ 
nance developed about four decades ago. The generalized 
cumulant method has proven to be a powerful tool for find¬ 
ing ensemble-averaged FID or solid echo components in the 
case when motions of the spin-bearing particles are uncorre¬ 


lated. In addition, the expansion of the density matrix in the 
set of symmetry-adapted eigenoperators allows one to take 
full advantage of the inherent symmetries of the problem, 
such as those associated with the statistical independence of 
the motions of spins and with the nonselectivity of pulses for 
like spins. 

We believe that the method of representing the commu¬ 
tation superoperator in an eigenoperator space may be ex¬ 
tended to more complex problems, including the relaxation 
of unlike spins and more involved pulse sequences. Among 
the advantages of this method are: more transparent symme¬ 
try properties of the interaction matrix and relatively simple 
vector-matrix multiplications involved. Moreover, it may al¬ 
low one to treat more complicated forms of the spin Hamil¬ 
tonian, such as for spins greater that 1/2. In subsequent work 
we intend to further explore the formal properties of the 
many-body interaction Hamiltonian in the eigenoperator rep¬ 
resentation for the purpose of developing a methodology for 
treating multiple quantum coherences and cases of unlike 
spins. 

The results in the present paper can best be studied in 
dilute solutions of like spins, of 1/2, where the dipolar inter¬ 
action is significant, and the slow motional regime may be 
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FIG. 6. Comparison of homogeneous and total broadening (. R 2 and R* , respectively) for multiple spins as a function of the diffusion coefficient D T , plotted 
on a logarithmic scale showing the solid state limit, the motional narrowing limit, and the intermediate motional region. Dashed lines indicate R* calculated 
using Eqs. (3.12) and (3.15) for the Redfield (fast motional) limit and the Anderson (solid-state) limit, respectively. The R 2 maximum occurs in the 
intermediate motional region; whereas in the motional narrowing limit R 2 and R 2 approach each other. 


conveniently achieved. This would be the case for ESR stud¬ 
ies on radicals, where slow motional effects set in even for 
moderately viscous fluids. Given the weak power depen¬ 
dence of Dj 032 for the homogeneous T 2 in the viscous (or 
slow motional) regime, this contribution would persist over a 
wide range of viscosities. An additional matter for study is 
the practical case of a complex spin Hamiltonian that in¬ 
cludes hyperfine couplings of the electron spin on each radi¬ 
cal to the magnetic nuclei, when present. This will require 
generalization of our approach to more complex spin Hamil¬ 
tonians along the lines suggested in the previous paragraph. 
Finally we note that only intermolecular spin interactions 
including intermolecular dipolar and (Heisenberg) spin ex¬ 
change terms require a many-body approach; intramolecular 
spin interactions, that are usually modulated just by rota¬ 
tional tumbling, can be more simply treated by the effective 
two-body stochastic motional analyses.of Paper I and refer¬ 
ences therein. 
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APPENDIX A: EIGENVALUES OF THE MATRIX SUMS, 

2" y [Cp]" 

The first important observation we shall make is that the 
right-hand side of Eq. (3.6) for g j(t) is invariant under a 
permutation of any two particles. On the other hand, permut¬ 
ing the two particles will at most switch two components of 
the left-hand sub-vector of coefficients g j(t) that have the 
same indices k and m. Making as many pairwise permuta¬ 
tions as necessary, one can then show that g[ k ' m \t) 
=g ( e k f ' m \t) for all e and e . In fact this is independent of 
whatever numbering system that is chosen to number the 
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permutation sets {e}. Therefore, it is sufficient to calculate 
only a single component of g/(/), i.e., g^ k ' m \t). 

Let us now consider the Ith block of the matrix 
defined in Eq. (2.10) and compute the result of the action of 
the matrix on the initial magnetization sub-vector, g/(0). For 
a given spin pair (ij) one can schematically represent the 
individual blocks of the matrix according to the classi¬ 
fication scheme for the coefficients g{*j w ^(0> cf. discussion 
below Eq. (2.16) noting that only coefficients g^' m \t) and 
g(*+i,m±i)(j) are C0U pj e( j by c] lJ \ Indices k and m will 
then vary from 0 to N— 1 —l and from / to 21—N+ 1, re¬ 
spectively. One can represent the product C/ ,j) g/(0) in the 
following form: 

(k, m) 



(Al) 

Here the nonzero elements of the matrix C\ 1 '^ are contained 
within the squares. Clearly, the result of the multiplication 
with blocks having k^2 will yield zeroes identically. There¬ 
fore, it is just necessary to discuss the blocks with k— 0 
and 1. 

In Eq. (Al) the diagonal block with k=0 will reflect 
couplings amongst operators of the same type, 
I + I a I a '-'l a IpIp"Ip, which contain exactly one operator 
/+, m = l operators I a , and N— 1— / operators Ip corre¬ 
sponding to the different spins. As follows from the commu¬ 
tation relations, cf. Eqs. (2.11), the block with k= 0 contains 
2 (N— 1) nonzero rows, which correspond to twice the num¬ 
ber of distinct pairs that one can form with the operators I + 
and the remaining (N— 1) operators I a and Ip. For a given 
spin pair (z'j), each nonzero row of the block will cortain 
exactly two numbers, ±1/2 and ±1/4, where the plus sign 
corresponds to the operation on I+^I^ or and a minus 

sign to the operation on I*+I*jp or I^I^p . Let us denote the 
positions of these two nonzero numbers (e,e) and re¬ 

spectively, where e is the index of the original coefficient 
g^' m \t) and e is the index of the coefficient g^,' m \t) that is 
obtained by switching of particles i and j. Now, the lower left 
off-diagonal (i.e., the [(0,/),(l,/)]) block in Eq. (Al) will 
reflect couplings between operators of two different types, 
namely, /+/+/_/*■ • -IJp- • 'Ip and I + IJ a ---I a IpIp---Ip. 
The operators of the former type having k= 1 can be num¬ 
bered by its own prescription {<5}. According to the commu¬ 
tation relations, each nonzero row of the off-diagonal block, 
e.g., rows £and S', will contain one +1/4 term and one —1/4 


term which will cancel each other upon multiplying by a 
column vector of all ones. Therefore, the result of action of 
the matrix on the initial magnetization vector g/(0) is, 



(A2) 

In order to compute [C^ 7y) ] 2 g/(0), we note that the col¬ 
umns of the lower off-diagonal block containing nonzero 
elements will always correspond to zero rows of the diagonal 
block of C^ /,7) with k = 0 [denoted as d' and d" in Eq. (A2)]. 
Indeed, as follows from the commutation relations, operators 
having k— 1 are obtained from operators with k — 0 when the 
combinations I^I^p ° r tyif are encountered, which cannot 
arise in the upper diagonal block describing the couplings 
between operators containing combinations I*+I^ and 
I^+I^jP [or I+^I^ and Z+V^]. Therefore, when multiplied by 
the matrix for the second time, the lower half of the 

resulting vector corresponding to £A0 will still consist of 
zeroes, since nonzero elements in rows 8 or S' of the lower 
off-diagonal blocks will be multiplied by zeroes at positions 
d' and d" arising from the previous action of the matrix 
C^. This means that the spin-flip terms in the dipolar 
Hamiltonian like I^I^P , apart from a numerical factor of 3/4 
in the interaction constant, do not affect the observable FID 
signal after a single nonselective tt/2 pulse. By contrast, the 
6th (or 6 r th) component of g/(0) with k— 0 will become 
(±3/4) 2 . Repeating this operation n times and summing 
over all spin pairs containing^ = / pairwise combinations 
I*+j^P each yielding one (3/4)", and ( N — 1 —m) combina¬ 
tions yielding one 

( — 3/4) w , one obtains for any 6 that 

{e [c^Tg/wj =[/(ir+(jv-i-/)(-!)”]. 

(A3) 

Applying the permutation symmetry argument given in the 
first paragraph of this Appendix, we conclude that g/(0) is an 
eigenvector of the matrix sums S / <y[C/ / ’- /) ]' 7 , with an eigen¬ 
value given by Eq. (A3). 
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APPENDIX B: EIGENVALUES OF THE MATRIX SUMS, 

2 IU C /'’ ;0 ] P X[Ci' v, ] n “ p 

We shall now prove that for a nonselective or hard pulse, 
each sub-vector of the starting vector g/(0) is an eigenvector 
of the matrix sum where / denotes 


the number of operators I a in an eigenoperator subset E^ 
having zero pairs 1+1- , i.e., &=0. 

As a result of a single nonselective — tt/2 pulse around 
the y axis applied at time r, the eigenoperators corre¬ 

sponding to fji= 1 will transform in the rotating frame as 
follows: 


L + a p + — 


- ](-ir/2) y 


2~ N (-I^ l) + I^ l) + / ( ± - I W ) 

X(^ 2) + / j2 ) + ^2 ) + /y)...(/03 ) + / U3 ) _^ 3) _ / (i3)). .. 


(- / ( ' 4) + / ( ‘ 4> +/ ( ± - / ( ' 4) ) ( - / (, ' s) + I { ‘ s) - 7 ( ^ 5) +/ (is) ) • • •, 

v a p + — a p + — / •> 


2k terms 


(Bl) 


r 


with an analogous expression for eigenoperators having 
jjl— — 1. Note that the signs of both the transformation to the 
rotating frame and the rotations by the pulse have been re¬ 
versed here to get directly the elements of matrix X, and not 
of its transpose, cf. Eq. (5.3). As can be seen from Eq. (Bl), 
a (tt72) v pulse mixes, in principle, all the components of the 
vector-coefficient g(t) with different coherence indices, jjl. 
However, in the high-temperature approximation, cf. Eq. 
(3.3), only the elements of the starting vector g(0) corre¬ 
sponding to k = 0 are nonzero, and, for the initial magnetiza¬ 
tion after the first 7t/2 pulse along the x axis, are equal (apart 
from the factor of2~ N h£ll2i JcT) to either +1 for fi= + 1, or 
— 1 for jul= — \, which significantly simplifies the calculation. 
According to Eq. (Bl), individual products containing spins / 
and j will yield the following combinations after the second 
pulse: 


r E 7 ') jUh 
l 1 + 1 a J(-tt/2) v 






’+/('■)/(/) 
— ±a 

-jU)M 

1 ± 1 a 

and 

and iflf. 

(B2a) 

-+T(i)jU) 
— ± 1 a 

i ±< ;) 

and ±I t ~+I^ 
and -iflp*. 

(B2b) 

rO) r U) 

-/</>/« 

and l ( i ] ] { p 
and -flf/ 

(B2c) 


and ±I^)l^p 
and + f+fp. 

(B2d) 


Other combinations containing l^I ^}, , and Ip*l*jp 

will yield zeroes upon further multiplication by the matrices 
[C^Y as follows from the commutation relations, Eq. 
(2.11), and, therefore, can be disregarded. The result of the 
matrix multiplication of the initial vector g(0) by matrices 
~ p is discussed in detail in Appendix A. In. particu¬ 
lar, for given k and m the vector g(0) contains the 


same values at positions corresponding to the permutations 
of the two spins (ij) and (/,/), i.e., (±3/4 ) n ~ p . When the 
vector [C^ J ^] n ~ p g(0) is further multiplied by X, for the 
positions in the resulting vector that correspond to combina¬ 
tions containing 1^+1 ^, the matrix multiplication over coef¬ 
ficients at 1^+1^ and 1+^1^ will yield zeroes; whereas the 
sums over positions corresponding to I^+l^jp or I+)lp* of 
[QOJ)y- p g(Q) w iH b e doubled. By contrast, for positions of 
the resulting vector pertaining to , the sums over po¬ 

sitions having 1+*I^ or I+^Ip* will vanish, while the sums 
over pairs I^+I^p and 1+^1^ will double. As can also be seen 
from Eqs. (B2c) and (B2d), the matrix multiplication over 
positions corresponding to or of the vector 

[C 0’J)y~p g(o) will always yield zeros. Therefore, for a 
given spin pair (ij) the result of action of the matrix X on 
vector g(0) will be just an exchange of values 

( — 3/4 ) n ~ p and (+ 3/4) n ~ p at the nonzero positions of the 
vector g(0). These values need to be further 

weighted by twice the sum over all permutations for the 
remaining N— 2 spins, including the sign of each permuta¬ 
tion with which it appears in Eq. (Bl). But this sum is for¬ 
mally just the sum over the coefficients of the polynomial 
represented in Eq. (Bl), which can be readily obtained from 
Eq. (Bl) by letting all the remaining N— 2 polarization op¬ 
erators l a and 1^ equal to unity and the remaining lowering 
and raising operators 1+ and /_ to be all zeros. Thus, for k 
= 0 this sum is equal to 2X2~ N X2 N ~ 2 = 1/2 and is zero for 
ki=- 0. Summing over the contributions from the components 
corresponding to both jx= + 1, having an additional factor of 
exp(—zAflr), and to jx= — 1 with exp(/Aflr), and invoking 
the permutation symmetry arguments similar to those used 
earlier to find eigenvalues of the matrix sum 
cf. Appendix A, one has for the 7th block of the coefficients 
gf^\t) that 
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N 

2 exp( — / AftT)&(0) 

><j 

= cos(A ClT)[ltin-\y-P 

+(7v-i-/)(-iFar^] g/ (o). (B3) 

As can be seen from Eq. (B3), it includes a factor of 
cos(AEIt) which diminishes the amplitude of the echo signal 
detected at the first-order coherence, since the remaining part 
of the signal is distributed amongst higher-order coherences. 
In this case, full forms of the matrices X, and All 

must be considered in order to derive Eq. (B3) rigorously, 
which will be a subject of the next paper. 
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